Impact of Exchange-Correlation Effects on the IV Characteristics of a Molecular Junction 
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The role of exchange-correlation effects in non-equilibrium quantum transport through molecular junctions is 
assessed by analyzing the IV curve of a generic two-level model using self-consistent many-body perturbation 
theory (second Born and GW approximations) on the Keldysh contour. For weak molecule-lead coupling we 
identify a mechanism which can lead to anomalously strong peaks in the dl / dV due to a bias-induced interplay 
between the position of the HOMO and LUMO levels. The effect is suppressed by self-interaction errors and is 
therefore unlikely to be observed in standard transport calculations based on density functional theory. Inclusion 
of dynamic correlations lead to substantial renormalization of the energy levels. In particular, we find a strong 
enhancement of quasi-particle (QP) scattering at finite bias which reduces the QP lifetimes significantly with a 
large impact on the IV curve. 

PACS numbers: 72.10.-d,71.10.-w,73.63.-b 
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Over the last decade it has become possible to contact sin- 
gle molecules by metallic electrodes and measure the IV 
characteristic of the resulting nano junction HI 0, 0] . Exper- 
iments of this kind can be seen as the first step towards the 
realization of a molecule based electronics. On a more fun- 
damental level, the IV characteristics provide a spectroscopic 
fingerprint of the molecular junction containing information 
about the positions and lifetimes of the electronic energy lev- 
els. In view of this, the interpretation of IV curves in terms 
of the electronic structure of the junction represents a funda- 
mental challenge for molecular electronics. 

So far, almost all ab intio calculations of conductance in 
molecular junctions have been based on the single-particle 
Kohn-Sham (KS) scheme of Density Functional Theory 
(DFT) El H. 

This approach has been successfully ap- 
plied to junctions characterized by strong coupling between 
molecule and leads |0 0], but has generally been found to 
overestimate the low-bias conductance of larger and more 
weakly coupled molecules 

HE EH. 

Recently, this short- 
coming of the DFT approach has been attributed to the pres- 
ence of self-interaction (SI) errors in the standard exchange- 
correlation (xc) functionals [l3ll . Inclusion of electronic 
correlations beyond the single-particle approximation could 
also be important, however, attempts in this direction have so 
far been limited 



EMU 



m\VM. 



While the theory of nano-scale conductance in the low bias 
limit and its relation to the ground state electronic structure 
has been studied in depth, the finite bias regime has received 



much less attention 11141 118L 11911 . The main reason for this 



unbalance is presumably the larger complexity of the latter 
problem: linear-response properties can be obtained from the 
ground state, e.g. via the Kubo formula [20], whereas finite- 
bias properties requires the construction of a non-equilibrium 
steady state which, unlike the ground state, is not a variational 
quantity. 

In this paper, we analyze the mechanisms governing the 
shape of molecular IV characteristics with particular focus on 
the role of exchange and correlation. The molecular junction 



is modeled by two electronic states representing the HOMO 
and LUMO levels symmetrically coupled to leads as sketched 
in Fig. [Tf a )- It is easv to anticipate that the slope of the IV 
curve is largest when a molecular level is aligned with one of 
the bias window edges. This will show as peaks in the dl/dV 
curve. In the simplest picture the distances between peaks in 
the dl/dV curve are thus a direct measure of the distances 
between the energy levels of the molecule, and the width of 
the peaks gives the lifetime of the levels. This simple picture 
breaks down because the electronic structure of the molecule 
will change in response to the applied bias, and as we will 
see, this effect can be surprisingly large even for very simple 
systems. 

On the basis of the two-level model we identify a simple 
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FIG. I: (a) Density of states for a molecular junction under zero and 
non-zero bias voltage. As indicated the bias changes the line shape 
of the DOS which in turn affects the dl/dV. (b) Diagrams for the 
second Born (2B) and GW self-energies. Full lines represent the 
Green's function of the molecule with coupling to leads. Wiggly 
lines represent the interaction. 
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mechanism which can lead to the formation of anomalously 
strong peaks in the dl/dV curve. The mechanism is driven 
by the applied bias voltage and can be viewed as a collapse 
of the HOMO-LUMO gap by which both levels move simul- 
taneously into the bias window giving rise to a large increase 
in the current. The mechanism is suppressed by SI errors and 
should therefore only be partly present in standard DFT trans- 
port calculations. The second result of the paper is that the 
inclusion of dynamic correlations, at the levels of the second 
Born (2B) and GW approximations, can have a strong influ- 
ence on the IV curve. This is mainly due to the enhancement 
of incoherent QP scattering at finite bias which leads a to sig- 
nificant broadening of the spectral features. This effect can 
also lead to a suppression of the anomalous dl/dV peaks. 

Our model consists of two electronic levels, representing 
the HOMO and LUMO states of the molecule, coupled sym- 
metrically to wide-band leads. The one-particle energies of 
the levels are £o and £o + Ao, where Ao is the non-interacting 
HOMO-LUMO gap. We take the charging energy of the 
HOMO and LUMO levels to be the same (Un) whereas the 
charging energy between an electron in the HOMO and an 
electron in the LUMO is set to U12 = 0.75Un. As discussed 
later our conclusions are not sensitive to this choice as long as 
U12 > 0.5t/ii. We neglect the exchange energy between the 
HOMO and LUMO states as it is generally much smaller than 
the charging energies [21]. The Hamiltonian of the molecule 
is written H mo \ = Ho + U, where 



= £]£o*W + (6) + Ao) 



(1) 



U = ^ Unh^hii + 22 U 12 nH a h La >, (2) 

i—H,L cr.er' 

with fiHa and Ulc being the number operator for an elec- 
tron with spin a in the HOMO and LUMO states, re- 
spectively. The interactions can also be written U = 

2 Ser,er' Yli je{L,H} Uia.jo'C^C^tCja'Cia, where Uicrj a i = 



Un5ij(l — <W') + ^12(1 — Sij). By working with this 
spin-dependent interaction, i.e. using U to represent the wig- 
gly lines in the diagrams of FigQJb), SI errors are automat- 
ically avoided to all orders in the interaction since Ui a ^ a — 
0. B22I1 This can be illustrated by noting that both the first- 
and second-order exchange diagrams (first and last diagram 
of T,2b) vanish when evaluated using ^7. 112311 In the present 
model these diagrams should exactly cancel the SI in the 
Hartree and second-order bubble diagrams, respectively. By 
using U for the calculation of the GW self-energy we get sim- 
ilar diagrams included "for free" to all orders in the interac- 
tion. In particular, the resulting Si-free GW is exact to second 
order, i.e. includes the 2B (at least for our model). This is 
in contrast to standard GW calculations, including previous 
GW model calculations 124 12511 . which suffer from higher 
order SI errors. j26ll 

The retarded Green's function describing the molecule in 
contact with leads is written 



Spin dependence has been suppressed as we specialize to the 
spin unpolarized case, i.e. Gjf = G^ — G l27tl . The cou- 
pling to leads is included via the wide-band tunneling rate T. 

and £ denote the Hartree and xc self-energies, respec- 
tively. In this work £ can be either exchange, 2B or (Si-free) 
GW. The self-energies, [G] and £[G], are calculated self- 
consistently in conjunction with the Dyson equation (0 using 
the non-equilibrium Keldysh formalism to account for the dif- 
ference in chemical potentials. The calculational procedure is 
described in detail in Ref . |22j. 

In the special case of symmetric coupling, the particle cur- 
rent can be written [28 ] 



I(V) = T / e {V;Lo)dto, 

Jep-V/2 



(4) 



G r {uj) = [{u + iT)I 



2x2 



H 



(3) 



where g(to) = (i/27r)Tr[G r (w) - G''(^) t ] is the (non- 
equilibrium) density of states (DOS). It has peaks at the posi- 
tion of the QP energy levels, which represent the electron 
addition/removal energies of the junction. The width of a peak 
equals the inverse lifetime of the QP, r^ 1 « r + ImE^). 
Because single-particle approximations are characterized by 
a real, frequency independent xc self-energies, the width of 
the spectral peaks obtained in Hartree and HF will be given 
by r. From Eq. © it is clear that the current will increase 
more rapidly when a peak in the DOS enters the bias window. 
We stress, however, that g{u) depends on V through £# and 
£, and thus the overall shape of the DOS, in particular the 
HOMO and LUMO positions, will change with the bias volt- 
age. Clearly, this change in the DOS determines the shape of 
the dl/dV curve. 

The following parameters have been used throughout: 
A = 2, Un = 2, U12 = 1.5. By varying the one-particle 
energy £ , we can control the equilibrium occupation of the 
molecule, N e \. We consider the case of weak charge trans- 
fer to the molecule, i.e. N e \ ranges from 2.0 to 2.1, corre- 
sponding to ep lying in the middle of the gap and slightly 
below the LUMO, respectively. The Fermi level is set to zero, 
and the bias is applied symmetrically, i.e. /jl = V/2 and 

m = -v/2. 

In Fig. [2] we show the calculated dl/dV curves (obtained 
by numerical differentiation) for different values of T and N e \. 
We first notice that the 2B and GW approximations yield 
similar results in all the cases indicating that the higher order 
terms in the GW self-energy are fairly small. For T = 1.0, all 
methods yield qualitatively the same result. For even larger 
values of T (not shown), and independently of N e \, the re- 
sults become even more similar. In this strong coupling limit, 
single-particle hybridization effects will dominate over the in- 
teractions. 

For all four sets of parameters, the Hartree approxima- 
tion severely overestimates the low-bias conductance. This 
is a consequence of the SI contained in the Hartree poten- 
tial which leads to an underestimation of the (equilibrium) 
HOMO-LUMO gap. On the other hand the HF, 2B, and GW 
methods lead to very similar conductances in the low-bias 
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FIG. 2: dl /dV curves for different values of the tunneling strength 
r and occupation of the molecule, 7V e i. The curves are calculated 
using different approximations for the xc self-energy. 



regime. This indicates that the inclusion of dynamic corre- 
lations does not change the (equilibrium) HOMO-LUMO gap 
significantly. It is well known that HF tends to overestimate 
band gaps, and that the inclusion of correlations, e.g. at the 
GW level, has the effect of reducing the HF gap towards the 
true value. However, in finite systems where the electrons are 
confined in discrete, well separated energy levels, correlation 
effects are weak and HF already yields good spectra. In such 
cases, the inclusion of correlations is expected to have minor 
effects on the equilibrium gap, as observed in the present case. 
We stress, however, that these considerations only apply to 
weakly coupled molecular junctions. In fact, for larger values 
of r the coupling to leads can provide significant screening of 
the interactions on the molecule, rendering the GW correction 
to the HOMO-LUMO significant J22]. 

Returning to Fig. [2] we notice that the lower left graph 
shows an interesting feature. Namely, the HF, 2B, and GW 
curves all contain an anomalously strong conductance peak. 
Interestingly, the peak height is significantly larger than 1 
which is the maximum conductance for a single level (the 
Anderson impurty model). Moreover, the full width at half 
maximum (FWHM) of the peak is only ohf = 0-27 and 
&2B/GW = 0.12, respectively, which is much smaller than 
the tunneling broadening of 2T = 0.5. We note in passing 
that the peak looses intesity as iV e i is increased, and that the 
Hartree approximation does not produce the anomalous peak 
at all. 



To understand the origin of the anomalous peak(s) we plot 
in Fig.[3]the evolution of the HOMO and LUMO positions as 
a function of the bias voltage (the 2B result is left out as it 
is similar to GW). Focusing on the upper panel of the figure 
(N e \ — 2.0), we notice a qualitative difference between the 
Hartree and the Si-free approximations: While the Hartree 



gap expands as the levels move into the bias window, the HF 
and GW gaps shrink leading to a dramatic increase in current 
around V = 2.5 and V = 1.3, respectively. This is clearly the 
origin of the anomalous dl/dV peaks. But why do the Si-free 
gaps collapse? 

Let us consider the change in the HOMO and LUMO po- 
sitions when V is increased by 25V. In general this change 
must be determined self-consistently, however, a "first itera- 
tion" estimate yields a change in the HOMO (LUMO) occu- 
pations of Sn H = -qh{-V/2)5V (Sn L = g L (V/2)SV). At 
the HF level this leads to 

de H = [-U lig (-V/2) + 2U 12 g(V/2)]5V (5) 
8e L = [U llQ {V/2) - 2U 12 g{-V/2)]SV (6) 

where we have used that qh{— V/2) ~ g(—V/2) and 
g L (V/2) w g{V/2). The factor 2 in front of U 12 includes 
interactions with both spin channels. In the symmetric case 
(N d = 2.0) we have q(-V/2) = g(V/2). Since U n < 2U 12 




FIG. 3: Position of the HOMO and LUMO levels as a function of 
the bias voltage for the Hartree (crosses), HF (triangles), and GW 
(circles) approximations. The horizontal lines show the FWHM of 
the GW resonances. The FWHM of the Hartree and HF resonances 
is 2T independently of V. Notice the differences in the way the levels 
enter the bias window: The Hartree gap opens while the HF and GW 
gaps close. In the upper graph T = 0.25, AT e i = 2.0 (symmetric 
case). In the lower graph V = 0.25, N e i = 2.1. 
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this means that 6eh > and Sel < 0, i.e. the gap is reduced 
as V is raised. Moreover it follows that the gap reduction is 
largest when g(±V/2) is largest, that is, just when the levels 
cross the bias window. In the general case (N e \ ^ 2.0) the 
direction of the shift depends on the relative magnitude of the 
DOS at the two bias window edges: a level will follow the 
edge of the bias window if the other level does not intersect 
the bias edge. It will move opposite to the bias, i.e. into the 
bias window, if the other level is close to the bias window 
edge. This is effect clearly seen in the lower graph of Fig. [5] 
(triangles). Thus the gap closing mechanism has the largest 
impact on the dl/dV curve when the HOMO and LUMO lev- 
els hit the bias window simultaneously. Moreover, the effect 
is stronger the larger Un/Uw, and the smaller T (the maxi- 
mum in the DOS is ~ 1 /T). At the Hartree level, Eqs. (0 and 
© are modified by replacing XJ\\ by 2XJ\\. This leads to an 
effective pinning of the levels to the bias window which tends 
to open the gap as V is increased, see Fig.|3](crosses). 

The effect of dynamic correlations can be identified by 
comparing the HF result with the GW/2B results in Fig. [2] 
For r = 0.25 two qualitative differences are observed: (i) The 
GW/2B conductance peaks occur at a lower bias voltage than 
the HF peak, (ii) The GW/2B peaks have long, flat tails on 
the high-bias side while HF peaks are more symmetric. Re- 
turning to Fig. |3]we see that part of (i) can be explained from 
the fact that HF yields a larger equilibrium gap than GW. In- 
deed, for V = the HF gap is - 0.3 larger than the GW 
gap. However, this effect alone cannot account for the shift 
in the conductance peak from V ~ 2.5 in HF to V ~ 1.5 in 
GW/2B. 

In fact, both (i) and (ii) are consequences of a significant 
spectral broadening occurring at finite bias in the GW/2B cal- 
culations. The broadening, indicated by horizontal lines in 
Fig. |3] is due to QP scattering. According to Fermi-liquid 
theory, QP scattering at the Fermi level is strongly suppressed 
in the ground state, i.e. Im£^(£j?) = for V = 0. However, 
as the bias is raised the phase space available for QP scatter- 
ing is enlarged and ImS increases accordingly. As a result 
of the additional level broadening, g(±V/2) increases more 
rapidly as a function of V. Since this is exactly the driving 
force behind the gap closing mechanism, the dl/dV peak oc- 
curs earlier in the GW than the HF calculation. Clearly, the 
long tails seen in the dl/dV of the GW 72B calculations are 
also a result of the spectral broadening due to QP scattering. 

In summary, we have investigated the mechanisms govern- 
ing the shape of the IV characteristic of a generic molecu- 
lar junction. We identified a simple gap-closing mechanism 
which can lead to anomalously strong peaks in the dl/dV. 
The mechanism is suppressed by SI errors and is therefore not 
likely to be correctly described in standard DFT calculations. 
This shows that the use of SI corrected xc-functionals, re- 
cently shown to be important for the linear-response conduc- 
tance, could be equally important under finite bias conditions. 
Finally, we found that the strong enhancement of QP scatter- 
ing as function of bias voltage leads to significant smearing 
of spectral features, which in turn has a large impact on the 



dl/dV curve. 
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